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ABSTRACT 


Very  often,  populations  exist  that,  logically,  should  satisfy  linear 
stochastic  ordering  requirements.  For  example  if  a  mechanical  device  is 
improved  through  N  stages,  the  corresponding  survival  functions  should 
he  linearly  stochastically  ordered.  Nevertheless,  estimates  may  not  re¬ 
flect  this  stochastic  ordering  because  of  the  inherent  variability  of  the 

observations . 

1.-  ■:  r'.v 

TTcrre  we-  characterize? the  maximum  likelihood  estimates  of  the  survival 

functions  subject  to  linear  stochastic  ordering  requirements.  —Wg~~sht>w -our 

estimates  may  be  expressed  in  terms  of  the  well-known  Kaplan-Meier  product 
. \  i  s 

limit  estimates.’  We-  also  give*, an  iterative  algorithm  which  we  show  must 

ft 

converge  to  the  correct  solution  that  depends  only  upon  solving  the  pair¬ 
wise  problem. 

G-it  !•  • 

Finally  we  consider  an  example  concerning  survival  times  for  people 
with  squamous  carcinoma  in  the  oropharynx  when  classified  by  degree  of 
lymph  node  deterioration  at  time  of  discovery. 
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1.  INTRODUCTION. 

Often  one  wishes  to  estimate  the  survival  functions  (l-cumulative 
distribution  function)  of  populations  from  possibly  censored  data  when 
nothing  is  known  concerning  the  functional  form  of  the  survival  functions. 
Appealing  estimates  in  this  situation  have  been  obtained  by  Kaplan  and  Meier 
(1958)  and  are  usually  referred  to  as  Kaplan-Meier  product  limit  estimators. 
Although  Kaplan  and  Meier  restricted  themselves  to  discrete  distributions, 
Johansen  (1978)  has  shown  that  the  product  limit  estimator  is  a  maximum  like¬ 
lihood  estimator  (mle)  in  the  class  of  all  distributions  under  the  generalized 
maximum  likelihood  framework  developed  by  Kiefer  and  Wolfowitz  (1956). 

Many  times  one  may  have  a  situation  where,  logically,  distributions 
must  be  stochastically  ordered.  For  example,  if  a  mechanical  device  is 
improved,  the  probability  of  survival  past  any  given  time 
for  the  improved  device  should  not  be  less  than  that  for  the  original 
device.  In  this  situation,  it  seems  reasonable  to  require  that  estimates 
of  the  survival  functions  should  also  satisfy  this  stochastic  ordering. 

Brunk  et  al.  (1966)  have  given  mle's  for  two  stochastically  ordered  cdf's 
for  uncensored  independent  random  samples.  Dykstra  (1982)  has  considered 
this  problem  for  the  case  of  right  censored  data  and  has  given  the  mle's 
in  the  form  of  Kaplan-Meier  product  limit  estimators  with  modified  data. 

In  this  paper,  we  are  able  to  find  the  mle's  of  NS  2  survival  functions 
when  a  linear  stochastic  ordering  exists  among  them.  While  the 
solution  has  a  nice  characterization, obtaining  the  actual  estimates  is 
quite  difficult.  An  iterative  algorithm  depending  only  upon  the  solution 
to  the  pairwise  problem  is  given  and  is  shown  to  converge  to  an  actual  mle. 


2.  THE  PROBLEM 


We  shall  assume  that  we  have  independent  random  samples,  possibly  with 
right  censored  observations,  from  N  discrete  populations.  The  N^2  popula¬ 
tions  are  assumed  to  be  stochastically  ordered,  so  that  the  corresponding 
survival  functions  satisfy 


(2.1) 


st  St 
P,  2  P  * 
1  2 


(We  say  P^  £  P^  if  P^(t)  i  P^(t)  for  all  t.)  The  problem  is  to  find 

nonparametric  mle's  of  the  survival  functions  subject  to  the  constraints  in 

(2.1).  We  will  without  loss  of  generality  (WLOG)  assume  that  P^O)  = 

P^O)  *  •**  =  =  1*  30  that  all  observations  are  positive. 

Complete  observations  (failures)  occur  on  a  subset  of  the  times 

S  <  S  <  •••  <  S  (S-  =Q,S  ,  =  “>).  The  number  of  failures  from  the  ith 
12  m  0  m+1 

population  which  occur  at  time  S.  is  denoted  by  d...  The  number  of 

J  1 J 

losses  (censored  observations)  in  the  interval  [S  ,S  )  from  the  ith 

J  J 

population  is  denoted  by  X  .  We  assume  the  X.  .  losses  occur  at  times 

ij  J 

,  r=l,*"*,/^  ,  where  these  censoring  times  are  fixed.  (The  same 

mle's  would  obtain  for  random  censoring  times  which  are  independent  of  the 

m 

times  of  failure.)  Let  n. .  =  Z  (d.  +X.  )  denote  the  number  of  items 

ij  r=tJ  ir  ir 

from  the  1  population  surviving  to  just  prior  to  S  .  We  have  assumed 
that  our  survival  functions  are  discrete,  but  this  is  not  really  necessary, 
as  one  may  argue  in  the  context  of  generalized  maximum  likelihood  (see 
Johansen  (1978))  that  our  estimates  need  place  probability  only  on  those 


timepoints  at  which  observations  occur. 


m 


! 
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3.  REDUCTION  OF  THE  PROBLEM. 

Based  on  the  notation  in  Section  2,  the  problem  is  to  find  survival 

functions  P  , • • • ,P„  which  maximize  the  likelihood 
1  N 


(3.1) 


N  fXiO 

IH  n 

i=llr=l 


p(L(i’0))  n  J 

2  j-ii 


[pi(sj-0)-pi(sj)] 


ij 


i.  . 

ij 

n 

r=l 


p  *  J  )  )l 

l  r  f 


subject  to  the  constraints 

(3.2)  P.(t)  i  Pi+1(t)  Vt,  i  = 1 » ' ' ' ,N-1 . 

For  a  given  set  of  survival  functions  P  ,***,Pjj  satisfying  (3.2), 
ve  note  that  the  likelihood  cannot  be  decreased  and  (3.2)  cannot  be  vio¬ 
lated  if  we  replace  P^(t)  by  a  discrete  P^(t)  which  has  possible 

Jumps  only  at  S, and  is  such  that  P^(S  )  =  P  (S  ) .  If  we  now 
t  m  1  J  1  3 

replace  P^(t)  by  P^(t)  defined  to  have  possible  Jumps  only  as 

S„  ,• • • ,S  and  P*(S.)  =  P-(S.),  the  likelihood  cannot  decrease  and  (3-2) 
1  m  2  j  2  J 

cannot  be  violated. 

Continuing  this  reasoning,  we  see  it  will  suffice  to  maximize  the 
expression 


N  m  d  X 

n  n  [P  (S  J  -P. (S  )1  P.(s  )  iJ 

i=l  J=1  J  J  J 

subject  to  P^S^)  a  Pi+l^Cj^  among  those  survival 

functions  which  place  probability  only  at  the  points  Si ,S2’ * " " ’^m‘ 


r 


h 


(We  note  that  maximum  likelihood  estimates  need  not  be  uniquely  defined 
if  the  last  observation  from  a  population  is  a  loss.  We  will  avoid  this 
ambiguity  by  requiring  our  maximum  likelihood  estimates  of  the  survival 
functions  to  be  as  small  as  possible.)  Equivalently,  we  wish  to  maximize 


N  m  r 

n  n  1 

i=l  j=lL 


Wi 


rP,<s4  J 


ij 


rr,  vd,  j  -iq  r  p.(s  )  -ix.. . 

.p^T  -  W]  [p^t  •••  Vsi>] 1J’ 


or  letting 


/ 


pi<s.i> 

TO' 


to  maximize 


N  m 

n  n  (i 

i=l  j=l 


n 

r<j 


pir 


i.  M.  , 
i.i  i.i 


subject  to  ft  p'  &  ft  p'  ,  j  =1»*  * "  »m;  i=l,***,N-l. 
r=l  lr  r=l  1  i’r 

m 

Finally,  recalling  that  n.  ,  *  Z  (d.  +£.  ),  letting  p.  .  =  In  v.  , , 

iJ  r_j  ir  ir  ij  ij 

and  considering  the  log  of  the  likelihood,  it  will  suffice  to  maximize 


(3.3) 


f(p 


1* 


N 

Z 


i=l 


m 

z 

J-l 


ln(l' 


pi  1 

■e  >*("u 


subject  to  the  constraints 


(3.U) 


j 

I  P 
r=l 


.  *  Z  pja,  ,  o  s  Pj  i 

ir  r=1  i+l»r  Fij 


for  j  =1, —  ,m;  i=l,***,N-l. 


The  problem  has  been  reduced  to  maximizing  a  concave  function  subject  to 
linear  inequality  constraints. 


5 


Since  we  are  maximizing  a  bounded  concave  function  over  a  closed 
convex  region,  there  must  exist  a  solution  p  =  (p^,***,p^)  which  satis¬ 
fies  the  constraints  in  equation  (3.^)  and  maximizes  equation  (3.3).  If 
the  active  constraints  (those  where  equality  holds)  of  equation  (3.M 
needed  for  a  solution  p  were  known,  the  number  of  variables  could  be 
reduced  by  using  certain  of  the  expressions  in  (3.^)  with  equality  holding. 
Then  expression  (3.3)  could  be  maximized  by  setting  the  appropriate  partial 
derivatives  equal  to  zero,  or  by  other  means,  and  solving  for  the  remaining 
independent  variables.  Of  course,  determining  which  constraints  are  the 
active  ones  is,  in  general,  a  very  difficult  problem.  However,  by  noting 
that  certain  equality  constraints  must  hold,  and  then  maximizing  equation 
(3.3)  subject  to  an  arbitrary  fixed  set  of  equality  constraints  in  (3.1*), 
the  general  form  of  the  solution  p  can  be  determined. 


a  a 

i  i 

a2 

a2 

Notice  that  the  constraints 

^  Pi2  ~  ^  Pi+1  2’ 

2=1  X  2=1  1,1 

al  ai 

E  Pi2  = 
2=1 

a2 

Ep. .=  Z  p. .  can  be  written  as  Z  p..=  Z  p. »,  Z  p. 
2= 1  lX  2=1  1+1  2=1  U  2=1  1+1 ’X  2=a,+l 


a2  \ 

^  ^i+1  l*  '  "*  ^  P-i  £  ~  *-• 

Z=a;L+l  1,X  X=a]c_l+1  1  X=a-  ’+1  1,1 


Z  p  so  that  an  arbitrary  p. 


need  be  present  in  at  most  two  active  constraints.  Thus  suppose  p. . 
is  present  in 


s  s— 1 

Pi,s  =  Pi-1,2  -  Pi  ,2 *  a  *  J  *  b  and  r  *  J  S  S- 

’  2=r  2=r 
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Also  suppose  that  cL^  >0,  and  no  other  parameter  equals  j •  Then, 

with  appropriate  substitution,  the  partial  derivative  with  respect  to 

p  of  equation  (3*3)  set  equal  to  zero  yields 
^  3 


(3.5) 


d  e  d 

- +  n,  -  d  =  h  -  k 

P, .  i 3  ij  i  j  i  j 

(l-e  lj) 


where 


,  Pi+l,b 

_  di+l,b6  ,  s 

'ij  ~  "  Pi+1  b  i+l,b"di+l,b) 

(l-e  1+1 


d.  e 

— — - +  ( n  -d .  ) . 

p  is  is 

(l-e  13 ) 


Of  course  it  doesn't  matter  which  variables  we  choose  as  the  depen¬ 
dent  variables.  It  follows  that  k. ,  must  have  constant  value  for 

ij 

a  s  j  s  b,  providing  d.  n  ,  >  0,  and  h. ,  must  have  constant  value 
for  r  «  J  s  s  if  d^j  >  0.  Solving  the  equation  for  p^  yields  a 
solution  of  the  form 


(3.6) 


/n, ,-d. ,+k. -h.  .  \ 

,*  .  JM  ij  xJ-iJ 

ij  V  nij+kij_hij  / 


ij  U  iJ 


Inspection  of  h^  and  k^  will  reveal  that  h^  =  k^_1  j  if  both 
d.  and  d.  n  are  positive.  If  d  =  0,  that  part  of  the  objective 

IJ  1-X, J  AJ 

function  (3.3)  involving  p^  is  linear  in  p^  with  nonnegative  slope. 
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Intuitively,  should  be  as  large  as  possible,  and  hence  equal  to  zero 

unless  our  constraints  will  not  allow  this.  This  would  happen  for  example 

if  d^^  >  0,  but  dg  =  d^g  =  0.  Then  our  constraints  would  require  that 

p.j^  =  p,^  and  p^^  =  P£2‘  ma^e  these  substitutions,  and  then  solve  by 

derivative  methods.  We  will  show  later  that  these  solutions  are  still  of 

the  form  (3.6)  for  appropriate  k.  ,  and  appropriate  definitions  of  inde- 

^  J 

terminate  forms.  Once  d. .  >0,  p*.  =  0  if  d.  .  =  0,  l  >  j.  Note  that 

ij  l  x*  ix 

if  p*  =  0  and  d.  =  0,  then  h.  ,  is  indeterminate,  and  hence  does 
Ut  iX  i  l 


not 


violate  our  claim  of  constant  h.  .  over  r  £  j  <.  s.  In  any  event,  p*. 

xj  xj 

can  always  be  written  in  the  form  of  (3-6)  if  indeterminate  values  are 
properly  defined. 

Since  the  true  solution,  p,  is  given  by  p*  for  some  set  of  active 
constraints,  p  is  necessarily  of  the  form 


PU  = 


In 


'nij~Vki.A 

nu+ku  ' 


(3.7) 


/n.  -d.  +k. .-k. 

/  ij  IJ  1.1  i- 

V  n. ,+k  -k.  , 

'  xJ  ij  x-1, 


p.  .  =  In 

ij 


=1A 

J 


Nj  \,  nNJ-kN_1^  J 


2  s  i  s  N-l 


for  appropriate  values  of  k^j  and  appropriate  definitions  of  indeterminate 
forms . 

Now,  the  problem  is  to  identify  the  k  's.  Heuristically ,  since  the 

-*■  J 

^ij's  may  be  interpreted  as  playing  the  same  role  as  the  n^'s,  the  esti¬ 
mates  are  obtained  by  transferring  unfailed  items  from  one  sample  to  the 
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next  stochastically  larger  sample,  and  then  using  the  Kaplan-Meier  (maximum 
likelihood)  estimate  for  each  survival  function.  It  will  then  turn  out 
that  the  k^'s  will  never  be  negative,  which  seems  reasonable  in  light 
of  our  one-sided  restrictions.  Of  course  probability  could  be  transferred 
from  one  population  to  several  larger  populations,  but  (3.?)  indicates 
that  we  do  not  need  to  think  of  things  this  way. 

Even  if  the  active  constraints  were  known,  a  large  number  of  simultaneous 
nonlinear  equations  would  have  to  be  solved  to  find  p.  In  the  following 
sections,  the  k^j's  will  be  characterized  and  a  method  using  this  char¬ 
acterization  will  be  given  to  estimate  the  N  survival  functions,  avoiding 
the  simultaneous  nonlinear  equations  problem. 
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k .  KUHN-TUCKER  CHARACTERIZATION 

The  characterization  of  the  k^j's  requires  some  theory  from  convex 
analysis.  The  problem  of  estimating  N  survival  functions  subject  to  the 
stochastic  orderings  given  by  equation  (3-2)  may  be  written  in  terms  of  a 
particular  ordinary  convex  programming  problem.  An  ordinary  convex 
programming  problem  is  to  find  the  vector  p  which  minimizes  a  given  con¬ 
vex  function  fp(p)  subject  to  the  constraints  f^(p)  s  0,  J  =  l,***,n  where 

R~ 

1 

our  problem  is  couched  in  the  terms  of  an  ordinary  convex  programming  prob¬ 
lem.  We  call  X  =  (X  ,X  , ***,X  )  a  Kuhn-Tucker  vector  if  X  £  0  for 

x  c  n  j 

j  =  l,***,n,  and  if  the  unrestricted  infimum  of  h(p)  =  fQ(p)  +  X^f^p)  + 

•••  +  X^f^p)  is  equal  to  the  restricted  infimum  of  f^(p).  The  following 
theorem  from  Rockafellar  (1970 )  characterizes  Kuhn  Tucker  vectors,  and  will 
be  useful  in  characterizing  the  k^^'s. 

Theorem  U.l.  Let  X  and  p  be  vectors  in  iP  and  ffi"1  respectively. 

~  „ 

In  order  that  X  is  a  Kuhn-Tucker  vector  for  our  problem  and  p  is  an 

optimal  solution,  the  following  conditions  are  necessary  and  sufficient. 

^  ^  ^  A 

a.  .X  £  0,  f  (p)  s  0,  and  X  f  (p)  =  0  for  j  =  l,---,n, 

J  J  J 

(U.l) 

b.  0  €  [SfQ(p)  +.  X^df^p)  +  •••  +.  X^Sf^p)] 

A 

where  dfj(p)  is  set  of>  subgradients  of  fj  at  p. 


each  f  is  a  convex  function  on  1 
J 


tion  (3.3),  and  letting  f  (p) 

i  j  j J=. 


-  i 


— - o 


1  •  ,m 


Since  the  true  solution  and  a  Kuhn-Tucker  vector  must  satisfy  the 
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subgradient  equation  in  (4. lb),  the  general  form  of  a  Kuhn-Tucker  vector 
can  be  found  for  our  problem  using  the  general  form  of  the  solution  given 
in  equation  (3.7). 


Theorem  4.2.  If  p  denotes  the  solution  to  the  problem  specified 
in  equations  (3-3)  and  (3-4),  then  a  Kuhn-Tucker  vector  corresponding  to 


the  constraints  f  *  E  p,  -  L  p. 

Jt=l  1  1,X  1=1 


-  ^  p,  ,  s  0  i 


s  given  by 


xu- 


Wj+i  •  1SJ<" 


when  the  solution  p  is  specified  in  the  form  given  by  (3.7). 


Proof;  In  the  ordinary  convex  program  characterization  of  our  problem, 

which  we  want  to  minimize  subject  to  the  constraints 

J  <3 

fii5(p)  *  I  pi+1  JL~  L  PiJt  s  0  for  i=l,***,N-l;  J  =  l,**,>m. 

jj=l  ji=l 

Notice  that  fQ(p)  and  f^j(p),  i  =  l,***,N-l;  j  =  l,»**,m  are  differenti¬ 
able  functions  of  p^  /  for  i*  and  j1  so  that  equation  (4.1b)  can  be 
written  in  matrix  notation.  (Of  course  the  n  in  Theorem  4.1  is  now 
m(N-l).)  Equation  (4.1b)  then  becomes  0  =  F_+FX,  where  F„  is  the 
NmX  1  vector  of  derivatives  of  fQ(p)  with  respect  to  p  , 


X  is  an 


(N-l)mXl  Kuhn-Tucker  vector  and  F  is  the  NmX  (N-l)m  matrix  of  par¬ 
tial  derivatives  of  f.  with  respect  to  p ,t  and  0  is  an  NmXl 

1  j  ~ 

vector  of  zeros. 

Equivalently,  we  must  have  Fk  =  -F  ,  and  since  F  is  of  full  column 

rank  which  is  less  than  the  number  of  rows,  any  solution  to  our  equation 

must  he  unique.  We  may  express  F  as 

0  0  •  •  •  0 

t  o  •••  : 

— T  •  •  • 

T  -T 
0  .  T 

where  T™011  is  an  upper  triangular  matrix  of  ones.  We  note  that  Fq  has 
the  form 


-k 


11 


-k 


lm 


k...  ~ 

ll  i+l,l 


k.  -  k 
im  i+l,m 


^-1,1 

kN-l,m 


where  i  =  1 ,  * • • ,N-2. 


It  can  he  verified  that  a  conditional  inverse  of  F  is  given  hy 
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Careful  calculation  will  verify  that  X  =  -FCF0  is  of  the  desired  form 
and  a  solution  to  equation  (4.1b). 

Conditions  (4.1a)  must  hold  if  the  X  's  are  to  form  a  Kuhn-Tucker 
vector.  If  X^  2:  0,  the  k^'s  occurring  in  equation  (3.7)  must  be 
such  that 

(4.2)  k.T  2  kJO  ^  ...  a  k.  2:  0  for  all  i. 

il  i2  xm 

Conversely,  if  equation  (4.2)  holds,  the  X.  will  be  nonnegative. 

-*-«J 

The  condition  that  X^f^(p)  =  0  in  Theorem  4.1,  together  with  the 
monotonicity  requirement  given  in  equation  (4.2)  leads  to  a  characterization 
for  the  solution  provided  certain  indeterminate  forms  are  properly  defined. 
Thus  we  have  the  following  theorem. 
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Theorem  4.3.  A  solution  to  the  problem  of  maximizing  f(p)  given  in 
equation  (3-3)  under  the  constraints  in  equation  (3.^)  are  given  by  the 
p^'s  described  in  equation  (3-7)  iff  the  k^'s  are  such  that: 

a)  the  p. ,'s  in  (3.7)  satisfy  the  constraints  in  (3.^)  and 
-‘s5ijs0- 

b)  kil  ^  ki2  s  ‘  £  kim  *  0  for  i  =  1,2,- • •  ,N-1 ,  and 


:)  “henever  = 


provided  some  indeterminate  forms  are  appropriately  defined. 


Indeterminate  forms  occur  when  d  =  0  for  iQ  <  i  and  J  s  and 

d.  >0.  In  this  case, 

0’J0 


PiQ, jQ  Pi0+l,J0  ’  PN,J0  ln 


Z  n  -d.  -k. 

l!=i0  1,J0  0,J0  1Q~1,JQ 

N 

Z  n.  ,  -k.  ,  . 

i=i0  1,J0  10~1,J0 


as  mentioned  earlier.  Notice  that  this  is  consistent  with  (3.7)  where  we 


kij  ”  i=J+1n*j  ’  1  >  i0’'5  S  j0 


since  this  leads  to  indeterminate  forms . 


5-  OBTAINING  THE  p 


Even  though  Theorem  U.3  characterizes  the  solution  to  the  problem, 
it  does  not  provide  any  means  of  obtaining  the  solution  in  general.  How¬ 
ever,  for  the  case  N  =  2,  Dykstra  (1982)  has  found  closed  form  expressions 
for  the  p  as  well  as  an  efficient  algorithm  for  calculating  them.  We 
will  elaborate  on  this  procedure,  and  then  give  an  iterative  procedure  based 
only  on  this  pairwise  scheme  which  will  converge  to  the  general  solution. 
This  will  give  us  a  means  of  solving  the  general  problem. 

To  solve  the  pairwise  problem,  one  need  first  identify  the  integer  JQ 
such  that  d  =  0,  j  <  j  ,  but  d  >0.  Then 


(5.1) 


plj  *  p2j  ■  ln 


ni.i  *ng.i  -ar 
"lj  +n2j 


(which  requires  that  )  for  j  <  jQ.  For  s  a  £  b,  we  let 

k  ,  denote  the  solution  to 
a,b 


b  [nn  -dn  +kl  b  /n?i  _d?  1  "k 

Z  IrJ-iJ - y—  -  z  - ^r— 

.  \  n..  ,  +k  I  ,  n0.  -k 

j=a  \  lj  /  j=a  2j 


if  it  exists  positively  and  0  otherwise.  Then  Dykstra  has  shown  the  p. 

J  j 

which  maximize  (3.3)  subject  to  the  linear  constraints  Lp,,  ^  L  p_ 

i=l  i=l 

are  of  the  form 


(5-2) 


pu  =  ln 


n„ ,  -d,  .  +  k. 


“u  +k« 


1  P2j  =  1 


n2j  "klj 


where  k,  .  *  min  max  k  .  for  j  a  j„.  The  restricted  mle’s  of  P.  (t) 
IJ  as  1  isb  a»b  0  x 


a£J  JSb  a»° 

and  P^(t)  are  then  given  by 
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(5.3) 


P.  (t)  -  exp  I"  E  p,  .1  , 

Lj;SjSt  1JJ 


i  =  1,2. 


We  shall  take  P^(t)  to  he  as  small  as  possible  subject  to  the  order 
restrictions  so  as  to  ensure  uniqueness  of  the  estimates.  Using  the  fact 
that  ^  a  •  •  •  k^,  Dykstra  has  derived  an  algorithm  which  can  be 

used  to  compute  the  mle  of  p  and  p  for  j  =  l,***,m.  This  algorithm 

J  *-«J 

is  given  below: 


Pairwise  Algorithm 


1. 

2. 


Algorithm  for  two  survival  functions  P^(t)  &  P^Ct): 

Define  p..^  (and  k  )  as  in  (5.l)  for  j  <  jQ. 

Find  the  largest  i  jQ  such  that  the  k^j  >  0  which  solves 

J1  /n  -d  +k\  J1  / n?,  -d  -k\ 

2  InMJ - ^r—  =  E  lnM^ - 1L- 

j=jQ  \  nlj  k  /  j=jQ  \  n2j  "k 


is  maximized.  Then  let  k., ,  =  k-  ,  for  jns  j  s  j  . 

JL  J  n 


*1J 


'0 


Find  the  largest  j  >  j  such  that  the  k  >0  which  solves 


In 


°i.t  -di.i  +* 


j-jj+i 


nii ♦* 


lnl!^-k' 


...  -  i  n_  ,  ■*  k 

j^j-j+l  \  2j 


is  maximized.  Let  k^  =  for  <  J  jg* 


4.  Continue  until  either  j  =  m  or  no  such  positive  k  .  exists,  in 


which  case  let  k,  =  0  for  j  >  j,  Then  p. . 

lj  h-l  ij 

given  as  in  (5.2)  and  (5-3)  respectively. 


lj 


and 


Pi(t)  are 
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The  N  sample  problem  can  be  solved  using  an  iterative  procedure 
based  upon  the  pairwise  algorithm  given  previously.  The  pairwise  method 
is  used  to  massage  (readjust)  the  data  until  an  equilibrium  point  is 
reached.  Letting  be  the  estimate  of  k  found  during  the  £— 

iteration  of  the  algorithm,  the  N  sample  algorithm  is  given  as  follows: 


1.  Find  the  k^j's  corresponding  to  P  s  P  using  the  pairwise  algo- 
rithm  given  previously. 

2.  Incorporate  kj1^  into  n  and  n  by  replacing  n  .  by 

““  n2,J  b> 

3.  Using  the  adjusted  data,  find  the  k^l's  corresponding  to  P  i  P 
using  the  pairwise  algorithm. 

Replace  n^-k^j  by  n2,  j +  “  kl  J  ”3,J  ^ 


5.  Continue  down  the  chain  of  survival  functions  until  the  N-l 
„th 


th 


and 


N  samples  have  been  compared,  and  n^  ^  has  been  replaced  by 
nK-l,J-i4-2,J*kH-l,J  “d  "»,j  »«  been  replaced  by  -  k'd>tJ . 


6.  The  procedure  starts  again  at  the  top  of  the  chain.  After  setting 
k^,  =  0  for  j  =  l,"’,m,  find  the  kj2|  corresponding  to  P  *  F 
using  the  pairwise  algorithm. 

(2)  (2) 

7.  Incorporate  k1  '  into  the  data  by  replacing  n  by  n  +k 


.(l) 

‘2,  j“2,j 


and  n^  ,+k^l  by  n0  , +  ki'"(-k 


U)  k(2) 

2,jT"2,j"kl,J- 


.(1)  . 


8.  After  setting  k'  '  =  0  for  J  =  l,***,m,  use  the  pairwise  algorithm 

2 ,  j 

(2)  (2) 

to  find  kl  ,  for  P  *  P,  and  incorporate  k0  ,  by  replacing 


k. 
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(2)  (2)  (2)  (l) 

n2,j-kl,j  by  n2J+k2,J-kl,j  “d  n3,j  +  k3,j  by 

n  +k(l)-k(2) 

3,j  k3,j  k2,j- 

9.  Continue  down  the  chain,  taking  out  the  effect  of  the  comparison 

(k!d^)  in  the  previous  iteration,  recomparing  P^  a  Pi+i>  an<i  'then 

(2) 

reincorporating  the  new  k^j  '  s. 

10.  When  the  bottom  of  the  chain  is  reached,  start  with  this  procedure 

again  at  the  top  of  the  chain  and  continue  until  achieving  convergence. 

The  estimated  p  's  are  then  obtained  from  (3.7)  and  the  latter  part 
^  J 

.  (jl) 

of  Theorem  4.3  with  the  k. ,  obtained  above  entered  in.  We  will  let 

ij 

indicate  (3.7)  with  the  kff^’s  entered, 
ij  U 


Theorem  3.1.  The  procedure  discussed  above  must  yield  values  of  p^ 
which  converge  for  every  i  and  j  as  4 Moreover,  these  limiting 
values  solve  the  problem  of  maximizing  f(p)  subject  to  the  constraints 

j  A 

Z  p.  .  *  I  p.A1  .  for  j  =  l,"*,m;  i  =  l,*",N-l. 

4=1  11  4=1  1+1 


Proof :  Consider  the  first  step  of  the  second  cycle.  We  obtain  k^. 

by  using  the  pairwise  algorithm  for  P^  i  P^  with  the  data  n ,  and 

n„.+ki'!')  in  place  of  n, ,  and  n„  .  We  need  to  consider  two  cases: 

2j  2j  r  lj  2j 

1.  If  d24=0  for  1SJ,  then  *  n2j*4j)  2  "2J  "  klj*' 

2.  Otherwise  k,  ,  =  min  max  k#  .  where  k*  .  is  positive  and 

lj  aij  jib  a»b  a>b 

solves 


d  /n  -  d  +  k\  d 
Z  In  MJ )  =  Z 

J=a  n  lj  /  j=i 


^  n2j  +  k2j  "k  ^ 
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or  is  zero. 


If  k&  ^  is  positive  and  solves 


,i  .1 


it  easily  follows  that  k7  ,  a  k  ,  for  all  a,b,  and  hence 

a,b  a,b 


k(2)  i  k(l) 
klj  ^  klj  * 


Similar  arguments  suffice  for  other  i,  and  we  may  conclude  that  k^  . 

(A) 

is  nondecreasing  in  A  for  all  i,j.  Since  the  k  are  uniformly 

N^-  J 

f  A) 

bounded  (say  by  E  n.J,  we  know  there  exists  k.  .  such  that  k.,  ^  k.  .. 

i=l  il  ij  i«5  iJ 

Finally,  noting  that  the  p  in  (3-7)  are  continuous  functions  of  the 

**■  J 

(A) 

k^  ,  and  that  the  k^  correspond  to  solutions  of  the  pairwise  problems, 
we  can  argue  that  the  p^  defined  by  (3.7)  for  the  limiting  values  k_^. 
satisfy  the  conditions  of  Theorem  U.3,  and  hence  solve  the  general  problem. 

We  mention  again  the  appealing  form  of  the  estimators.  They  are  still 
of  the  Kaplan-Meier  form,  where,  heuristically ,  unfailed  items  from  a 
population  are  transferred  to  the  next  stochastically  larger  population. 
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6.  EXAMPLE 

This  example  involves  a  large  clinical  trial  carried  out  by  the 

Radiation  Therapy  Oncology  Group  cited  in  Kalbfleisch  and  Prentice  (1980). 

The  analysis  was  done  on  only  a  small  part  of  the  data,  and  females  were 

deleted  to  make  the  data  more  homogeneous.  Patients  diagnosed  with 

squamous  carcinoma  of  the  oropharynx  were  classified  by  the  deg.ee  to 

which  the  regional  lymph  nodes  were  affected  by  this  disease. 

We  let  P  be  the  survival  function  corresponding  to  those  patients 

for  whom  the  disease  had  no  effect  on  the  lymph  nodes.  Correspondingly, 

we  let  P^,  P^,  and  P^  be  the  survival  functions  of  patients  with 

increased  effects  of  the  disease  on  the  lymph  nodes.  The  data  for  each 

survival  function  is  given  in  Table  1. 

Because  of  the  nature  of  the  disease,  one  would  certainly  expect  a 

linear  ordering  with  respect  to  lymphatic  involvement.  Thus  an  ordering 
st  st  st 

of  the  form  ^  P2  4  P3  ^  seems  inherently  reasonable. 

The  unrestricted  mle's  (Kaplan-Meier )  for  these  survival  functions 

are  given  in  Table  3.  Notice  that  with  respect  to  the  expected  orderings, 

there  are  reversals  so  that  estimates  requiring  these  orderings  will  be 

different.  The  algorithm  given  in  Section  5  was  applied  to  the  problem. 

The  k..  values  found  in  the  first  two  iterations,  in  the  order 

ij 

in  which  they  were  found,  along  with  the  essentially  limiting  values 
obtained  after  20  cycles  are  given  in  Table  2.  (The  procedure  used  here 
actually  worked  pairwise  from  bottom  to  top  rather  than  from  top  to  bottom 
as  indicated  in  the  paper.  Of  course,  the  limiting  values  will  be 
identical.)  The  values  of  these  kj^  indicate  how  much  "smoothing"  is 

*J 

occurring  between  populations. 


The  mle's  subject  to  the  stochastic  ordering  P^  i  i  P3  S  P^, 
obtained  via  our  algorithm  are  displayed  in  Table  U.  Note  that  the  sto¬ 
chastic  ordering  requirement  has  brought  substantial  change  in  our 
estimates. 
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TABLE  1 

Survival  Times  in  Days  of  Patients  with  Squamous  Carcinoma  in  the  Oropharynx 
with  Various  Degrees  of  Lymph  Node  Deterioration  (+  denotes  censored  observation). 


Degrees  of  deterioration  of  lymph  nodes. 


i=  1 

i  =  2 

i  =3 

i=  1* 

167 

216 

105 

9k 

238 

32l* 

222 

99 

276+ 

338 

279 

112 

296 

3U7 

395 

127 

32k 

599 

1+65 

131+ 

351 

763 

5  h6 

11+7 

372 

929 

915 

192 

37 k 

1086+ 

918 

219 

hob 

1092 

1058+ 

255 

1+1+5+ 

1317+ 

1U55+ 

262 

5kl 

1609+ 

161*U+ 

271+ 

560 

293 

9U3 

307 

998+ 

327 

123U+ 

33l+ 

ll+6o+ 

363 

1823+ 

1*07 

1+13+ 

(The  degree  of 

deterioration 

( i )  equals 

1+59 

517 

N-stage  tumor 

classification 

in 

532 

Kalbfleisch  and  Prentice  (1980)  plus  one.) 

51+1+ 

672 

696 

800 

9lU 

1312+ 

11+1+6+ 

11+72+ 
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TABLE  2 


Values  of  the  k. 


obtained  for  survival  data  in  Table  1. 


First  cycle: 


(l)  _  2.333 


j  ^  3 
j  >  3 


T  £  105 
T  >  105 


(1)  X-853 
krj;  =  1.205 

J  0 


3  *  29 
29  <  j  *  53 
j  >  53 


T  £  37*+ 
37 1+  <  T  £  1609 
T  >  1609 


(l)  12.853 

kUj=  U.380 
ij  0 


3  *  9 
9  <  3  £  1+0 
3  >  *+0 


T  £  192 
192  <  T  £  560 
T  >  560 


Second  cycle: 


(2)  _  3.569 
3J  0 


j  *  3 
j  >  3 


T  £  105 
T  >  105 


5.1+81 
(2)  _  3.752 
L2j  1.205 
0 


j  *  11 
11  <  j  £  29 

29  <  3  *  53 
j  >  53 


T  £  219 
219  <  T  £  371+ 
37l+  <  T  £  1609 
T  >  1609 


t9\  12.853 

k \2.>  =  5.8U7 
XJ  0 


j  *  9 
9  <  j  as  1+0 
j  >  1+0 


T  £  192 
192  <  T  £  560 
T  >  560 


Twentieth  cycle: 

(  .  10.525 

kV;0)  =  1.1+38 
3j  0 


3  *  3 
3  <  j  £  17 
J  >  IT 


T  £  105 
105  <  T  £  279 
T  >  279 


12.290 
(20)  _  5.763 

L2j  1.205 
0 


j  *  11 
11  <  J  £  29 
29  <  j  £  53 
j  >  53 


T  £  219 
219  <  T  £  37l+ 
37l+  <  T  £  1609 
T  >  1609 


12.853 

10.371 

7.253 

0 


j  *  9 

9  <  j  £  20 
20  <  j  £  1+0 

j  >  ^0 


T  £  192 
192  <  T  £  307 
307  <  T  £  560 
T  >  560 


TABLE  3  23 


Kaplan-Meier  MLE's  of  survival  functions  from  data 


Index 

Time 

Pl 

P2 

P3 

P4 

0 

0.0 

1.0000 

1.0000 

1.0000 

1.0000 

1 

94 

1.0000 

1.0000 

1.0000 

.9655 

2 

99 

1.0000 

1.0000 

1.0000 

.9310 

3 

105 

1.0000 

1.0000 

.9091 

.9310 

4 

112 

1.0000 

1.0000 

.9091 

.8966 

5 

127 

1.0000 

1.0000 

.9091 

.8621 

6 

134 

1.0000 

1.0000 

.9091 

.8276 

7 

147 

1.0000 

1.0000 

.9091 

.7931 

8 

167 

.9412 

1.0000 

.9091 

.7931 

9 

192 

.9412 

1.0000 

.9091 

.7586 

10 

216 

.9412 

.9091 

.9091 

.7586 

11 

219 

.9412 

.9091 

.9091 

.7241 

12 

222 

.9412 

.9091 

.8182 

.7241 

13 

238 

.8824 

.9091 

.8182 

.7241 

lU 

255 

.8824 

.9091 

.8182 

.6891 

15 

262 

.8824 

.9091 

.8182 

.6552 

16 

274 

.8824 

.9091 

.8182 

.6207 

IT 

279 

.8824 

.9091 

.7273 

.6207 

18 

293 

.8824 

.9091 

.7273 

.5862 

19 

296 

.8193 

.9091 

.7273 

.5862 

20 

307 

.8193 

.9091 

.7273 

•  5517 

21 

324 

.7563 

.8182 

.7273 

.5517 

22 

327 

.7563 

.8182 

.7273 

.5172 

23 

334 

.7563 

.8182 

.7273 

.4828 

24 

338 

.7563 

.7272 

.7273 

.4828 

25 

347 

.7563 

.6364 

.7273 

.4828 

26 

351 

.6933 

.  6364 

.7273 

.4828 

2T 

363 

.6933 

.6364 

.1213 

.4483 

28 

372 

.6303 

.6364 

.7273 

.4483 

29 

374 

.5672 

.6364 

.7273 

.4483 

30 

395 

.5672 

.6364 

.6364 

.4483 

31 

4o4 

.5042 

.6364 

.  6364 

.4483 

32 

407 

.5042 

.  6364 

.6364 

.4138 

33 

459 

.5042 

.6364 

.  6364 

.3762 

34 

465 

.5042 

.  6364 

.5455 

.3762 

35 

517 

.5042 

.  6364 

.5455 

.3386 

36 

532 

.5042 

.  6364 

.5455 

.3009 

37 

541 

.4322 

.6364 

.5455 

.3009 

38 

544 

.4322 

.6364 

.5455 

.2633 

39 

546 

.4322 

.  6364 

.4546 

.2633 

4o 

560 

.3601 

.6364 

.4546 

.2633 

hi 

599 

.3601 

,5455 

.4546 

.2633 

42 

672 

.3601 

.5455 

.4546 

.2257 

43 

696 

.3601 

.5455 

.  4546 

.1881 

kk 

763 

.3601 

.4546 

.4546 

.1881 

800 

.3601 

.4546 

.4546 

.1505 

U6 

914 

.3601 

.4546 

.4546 

.1129 

**T 

915 

.3601 

.4546 

.3636 

.1129 

»*8 

918 

.3601 

.4546 

.2727 

.1129 

49 

929 

.3601 

.3636 

.2121 

.1129 

50 

943 

.2881 

.3636 

.2121 

.1129 

51 

1092 

.2881 

.2424 

.2121 

.1129 

52 

1472 

.2881 

.2424 

.2121 

.0000 

53 

1609 

.2881 

.0000 

.2121 

.0000 

54 

1644 

.2881 

.0000 

.0000 

.0000 

55 

1823 

.0000 

.0000 

.0000 

.0000 

l 


MLE's  of  survival  functions 


Index 


TABLE  U 
under  order 
P, 


restrictions  from  data  in  Table  1. 


0.0 
9h 
99 
105 
112 
127 
13l 
lbl 
167 
192 
216 
219 
222 
238 
255 
262 
27k 
219 
293 
296 
307 
32h 
327 
33b 
338 
3^7 
351 
363 
372 
31 b 

395 

i*oi» 

1+07 

1*59 

1+65 

517 

532 

5^1 

5I+U 

5U6 

560 

599 

672 

696 

763 

800 

91b 

915 

918 

929 

9b3 

1092 

11*72 

1609 

161*1* 

1823 


1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

.91*59 

1.0000 

1.0000 

1.0000 

.8918 

1.0000 

1.0000 

.8917 

.8918 

1.0000 

1.0000 

.8917 

.8569 

1.0000 

1.0000 

.8917 

.8220 

1.0000 

1.0000 

.8917 

.7871 

1.0000 

1.0000 

.8917 

.7522 

.9665 

.9665 

.8917 

.7522 

.9665 

.9665 

.8917 

.7173 

.9665 

.8917 

.6917 

.7173 

.9665 

.8917 

.8917 

.6821* 

.9665 

.8917 

.731*6 

.6821+ 

.9299 

.8917 

.731*6 

.6821* 

.9299 

.8917 

.731+6 

.61*76 

.9299 

.8917 

.7316 

.6127 

.9299 

.8917 

.7316 

.5778 

.9299 

.8917 

.5775 

.5116 

.9299 

.8917 

.5115 

.  5U57 

.8917 

.8917 

.5115 

.51*57 

.8917 

.8917 

.5115 

.5136 

.81*77 

.7869 

.5115 

.5136 

.81*77 

.7869 

.5115 

.1*815 

.81*77 

.7869 

.5115 

.  1*1*9!+ 

.81*77 

.6821 

.5115 

.1*1*91* 

.81*77 

.5773 

.5115 

.1*1*91* 

.8036 

.5773 

.5115 

.1*1*91* 

.8036 

.5773 

.5115 

.1*173 

.1596 

.5773 

.5115 

.1*173 

.7156 

.5773 

.5115 

.1*173 

.7156 

.5773 

.1*925 

.1+173 

.6716 

.5773 

.1*925 

.1*173 

.6716 

.5773 

.1*925 

.3852 

.6716 

.5773 

.1*925 

.3502 

.6716 

.6716 

.6716 

.621*5 

.621*5 

.621*5 

.5113 

.5773 

.5773 

.5773 

.5773 

.5113 

.5773 

.5773 

.5773 

.5773 

.1*619 

.1*619 

.1*619 

.1*619 

.1*619 

.0000 


.5773 

.5113 

.5773 

.5773 

.5773 

.5773 

.5773 
.5070 
.5070 
.5070 
.1*366 
.  1*366 
.1*366 
.  1*366 
.  1*366 
.3662 
.3662 
.2791 
.2191 
.1525 
.0000 
.0000 


.1+075 

.1*075 

.1*075 

.1*075 

.1+075 

.3225 

.3225 

.3225 

.3225 

.3225 

.3225 

.3225 

.3225 

.2375 

.1526 

.1526 

.1526 

.1526 

.1526 

.1526 

.0000 

.0000 


.2801 

.2801 

.21*51 

.21*51 

.21*51 

.21*51 

.2101 

.1751 

.1751 

.11*01 

.1051 

.1051 

.1051 

.1051 

.1051 

.1051 

.0000 

.0000 

.0000 

.0000 
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Very  often,  populations  exist  that,  logically,  should  satisfy  linear  sto¬ 
chastic  ordering  requirements.  For  example  if  a  mechanical  device  is  improved 
through  N  stages,  the  corresponding  survival  functions  should  be  linearly 
stochastically  ordered.  Nevertheless,  estimates  may  not  reflect  this  stochas¬ 
tic  ordering  because  of  the  inherent  variability  of  the  observations. 

Here  we  characterize  the  maximum  likelihood  estimates  of  the  survival 
functions  subject  to  linear  stochastic  ordering  requirements.  We  show  our  esti' 
mates  may  be  expressed  in  terms  of  the  well-known  Kaplan-Meier  product  limit 
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estimates.  We  also  give  an  iterative  algorithm  which  we  show  must  converge 
to  the  correct  solution  that  depends  only  upon  solving  the  pairwise  problem. 

Finally  we  consider  an  example  concerning  survival  times  for  people 
with  squamous  carcinoma  in  the  oropharynx  when  classified  by  degree  of  lymph 
node  deterioration  at  time  of  discovery. 
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